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Abstract. We apply an ad hoc model for dynamical ignition in three-dimensional numerical simulations of thermonuclear 
supernovae assuming pure deflagr ations. The m odel makes use of the statistical description of temperature fluctuations in the 
pre-supernova core proposed bv lWunsch & WooslevI i2004) . Randomness in time is implemented by means of a Poisson pro- 
cess. We are able to vary the explosion energy and nucleosynthesis depending on the free parameter of the model which controls 
the rapidity of the ignition process. However, beyond a certain threshold, the strength of the explosion saturates and the out- 
come appears to be robust with respect to number of ignitions. In the most energetic explosions, we find about O.TSM© of iron 
group elements. Other than in simulations with simultaneous multi-spot ignition, the amount of unburned carbon and oxygen at 
radial velocities of a few 10^ km/s tends to be reduced for an ever increasing number of ignition events and, accordingly, more 
pronounced layering results. 
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1. Introduction 

The significance of the ignition characteristics at the on- 
set of thermonuclear supernova explosions has been recog- 
nised since the advent of multi-dimensional numerical simu- 
lations. The mounting evidence for diversity among observed 
type la supernovae renders the quest for realistic mod elling 
of the ignition proc e ss yet more topical (iBenetti et al.l |2005 ; 
lLeonard etal.1 120051: IStritzinger et all [2005"). Substantial ef- 
forts have been made to analyse the eflTect of choosing dif- 
ferent initia l burning regions in simulations of thermonuclear 
super novae ( Garcia- Senz & Bra vo '20051 iRopke & HillebrandJ 
120051) . The multi- point ignition scenarios which have recently 
been presented by iRopke etahl (l2005l) suggest that the energy 
generation becomes maximal for about one hundred ignitions. 
In this case, unburned carbon and oxygen at low radial veloc- 
ity largely diminis h. Except for the two-dimensional study by 
iLivne et all fe005'), however, it has not been attempted so far to 
generate ignitions at diflferent instants of time during a simula- 
tion. 

In this article, we present a simple stochastic model for 
the dynamical ignition of thermonuclear deflagration in three- 
dimensional simulations. The purpose of the model is an in- 
vestigation of the possible consequences of using a simple 
parametrisation for the ignition probability. In particular, the 
total amount of iron group and intermediate mass elements as 
well as their spatial distribution and density functions in radial 
velocity space are studied. We do not claim yet that the highly 
complex physics of the ignition process in thermonuclear su- 
pernova is fully taken into account. Notwithstanding the ten- 
tative nature of the stochastic ignition model, it may help to 



clarify the extent of possible variations in the yields and the 
distribution of burning products without invoking a hypothet- 
ical d eflagration to detonation transition (iGamezo etal .112004 
I2OO5I) . Whether the predictions of the model can be reconciled 
with properties of observed type la supernovae (SNe la) will 
depend on further developments in the theoretical understand- 
ing of convection and the birth of flames in the pre-supernova 
core. 



2. Stochastic ignition 

The simplest thermodynamical description of the pre- 
supernova core is an adiabatic background plus temperature 
fluctuations due to convection driven by gradual thermonu- 
clear burning. For the background temperature, we assume the 
parabolic profile 



ra(r) = max < T( 



(1) 



with a cutoflT at the thermal radius A = 7.35-10^ cm for the 
central density pc = 2 • 10^ g cm"^ (IWunsch & Wooslevll2004l) 
and an isothermal exteri or region with Text = 5 • 10^ K. For the 
temperature fluctuations. IWunsch & WooslevI (l2004l) proposed 
a statistical model on the basis of the mixing-length approach. 
The typical temperature excess dT\y gained by fluid elements 
relative to the adiabatic background in a spherically symmetric 
core region with Kolmogorov power spectrum is given by 



-3/7 
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where R can be interpreted as mixing length, is the loga- 
rithmic derivative of density relative to pressure a nd r is the 
nuclear time scale (see IWunsch & WooslevlEo04[ eqn. 9 & 
18). 

The ignition algorithm utilised in our simulations is based 
on a discrete stochastic process, namely, the Poisson process. 
The probability for / ignitions within the time interval St is 
given by 



(VtrSty 



(3) 



where Vt^r is a function of the temperature excess ST\) and the 
local thermodynamic background state. Assuming that the ig- 
nition events are statistically independent, the number of igni- 
tions per unit time increases proportional to the volume of fuel 
within a region of approximately constant background con- 
ditions. For this reason, we consider the ignition probability 
Pt,r,6r(i, St) for the unburned material within spherical shells of 
width Sr ^ R. On grounds of ergodicity, the increase of prob- 
ability in proportion to the volume is equivalent to either in- 
creasing the sampling time or decreasing Vt^r by the ratio 
of the fuel volume relative to a unit volume. Thus, Pt,r,6r(i, ^0 
is given by an expression analogous to the right hand side in 
equation (|3 with replaced by 



yt,r,6r = {yt,rO[G(t,r)])r, 



3r Sr 



.Sr- 



(4) 



where G(t,r) is the level set function which represents the 
flame front, 6[G(t, r)] is zero for G(t, r) < (i.e. in ash regions) 
and unity for G(t, r) > 0. ( )r,dr denotes the average over a 
spherical shell of radius r and width dr. For the normalisation, 
we choose a sphere of radius ry, which is roughly t he size of 
the central heating region ( W unsch & Woosle vl2004 eqn. 24). 
The underlying assumption of spherical symmetry is, of course, 
increasingly flawed in the course of the explosion. However, we 
found empirically that the spatial variation of the background 
state variables within the unburned material contained in a thin 
spherical shell at each given instant of time is small. 

In order to complete the model, we formulate a hypothesis 
for {vt,rO[G{t,r)\)r,dr- Since Vt,r is an inverse time scale which 
determines the probability of ignition per unit time, we conjec- 
ture that 

Xu,r,6r{tr\yt,re[G{U r)])r,6r = ^e^^, (5) 

where is the heat capacity and 5'nuc is the nuclear energy 
generation rate for the average temperature and mass density, 
respe ctively, of unburned m aterial within a particular spherical 
shell ("Woosley et al. 2004, eqn. 4 and 6). The exponentiation 
parameter Ce determines the overall rapidity of the ignition pro- 
cess. The filling factor ;^u^^^^^(0 specifies the volume fraction of 
unburned fluid at time t within the shell. Note that both sides 
of equation Q are independent of the amount of fuel which is 
consistent with the statistical nature of the quantities involved. 

3. Numerical simulations 

For the numerical simulations presented in the following, we 
use the same techniques as described in >Ropke & Hillebrandt 



(EOO?), except for a hybrid grid geometry with a central 
uniform part and exponentially increasing grid size outside 
(iRonke etZll2005h in combination with a localised subgrid 
scale model which is introduced in ISchmidt et alJ (l2005h . 
Owing to the co-moving hybrid grid, we were able to achieve 
sensible results with merely 128^ numerical cells in one oc- 
tant. The subgrid scale model is based on the hydrodynami- 
cal equation for the turbulence energy which was adopted by 
iNiemever & Hillebrandtl (Il995l) for the numerical simulation 
of thermonuclear supernovae. In the new implementation, we 
apply dynamical procedures for the in situ calculation of clo- 
sure parameters. Moreover, an additional term was added to the 
subgrid scale turbulence energy equation which accounts for 
buoyancy efl'ects on unresolved scales ( Schmidt et al. 2005). 
Thereby, the Sharp-Wheeler relation for the characteristic ve- 
locity scale of the Rayleigh-Taylor instability is efl'ectively in- 
corporated into the turbulence energy model. In addition, the 
localised SGS model features an algorithm which corrects the 
SGS turbulence energy for the efl'ect of the shifting cutoflf 
length scale of the co-moving grid, because the gradually in- 
creasing cutofl' length entails an intrinsic growth of SGS turbu- 
lence energy ( Schm idt et al. 2005). 

The nume rical im plementation of the Poisson process was 
adopted from Press et alJ 41296). For each ignition event in a 
particular Poisson process realisation, a random point is sam- 
pled from a uniform spatial distribution and projected onto the 
spherical shell associated with the process. If the resulting point 
happens to fall inside fuel, a small bubble of burned material 
with radius 1.5A, where A is the numerical cell size, is inserted. 
This algorithm implicitly accounts for the filling fdiCiov X\x,r,6rif) 
of the fuel in each shell. For the width of the spherical shells 
we chose Sr = 2A. 

We performed simulations with diff'erent values of Ce both 
in octants and full spherical domains. The influence of difl'erent 
mixing lengths R relative to the thermal radius A was investi- 
gated as well. The simulation parameters as well as the total 
number of ignitions, the explosion energetics and the burning 
product masses are listed in Table [T] In one series of simula- 
tions the value of the exponentiation parameter Ce was var- 
ied over four orders of magnitude. The main conclusion is 
that the explosion energy and the total mass of iron group el- 
ements (represented by ^^Ni) is maximal for Ce ~ 10^. If the 
ignition process becomes even more rapid, bubbles of burn- 
ing material are generated at such high rate that most of the 
bubbles merge before they can appreciably expand and deform 
due to the Rayleigh-Taylor instability. This efl'ect tends to in- 
hibit the growth of the flame surface area and, thus, reduces 
the explosion energy. More or less the same trend was ob- 
se rved for simultaneo us ignition with varying bubble number 
bylRopke etalJilooi. 

For Ce <^ 10^, there are only few ignitions per octant. In 
this case, the explosion is weak and the total mass of iron group 
elements is less than 0.5Mq. Setting Ce = 10, the white dwarf 
becomes barely unbound. Nevertheless, the ignition process 
ceases because the probability of further ignitions becomes 
virtually zero once the star is expanding appreciably. This is 
simply a consequence of the power-law dependence of the fac- 
tor Snuc/cp in formula ^ on the background temperature and 
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Table 1. Parameter study with initial spatial resolution A = 4.75 km in the uniform part of the numerical rid. The total spherical 
angle covered by the simulation domain is fl, the total number of ignitions per octant is /;r/2, and the last ignition event occurs at 
time ^lasf 





Ce 


R/A 


Q 


In/2 


hast [s] 


^nuc[10^^erg] 


£kin[10^^erg] 


MkJMq 


Mug/MQ 




10^ 


0.5 


7T/2 


2 


0.200 


0.516 


0.008 


0.242 


0.165 




10^ 


0.5 


7T/2 


4 


0.096 


0.827 


0.320 


0.413 


0.218 




103 


0.5 


tt/? 
71/ Z 


30 


0.508 


1.139 


0.631 


0.629 


0.187 


5- 


10^ 


0.5 


7T/2 


65 


0.177 


1.219 


0.711 


0.690 


0.170 




10^ 


0.25 


n/l 


65 


0.199 


1.099 


0.591 


0.603 


0.189 




10^ 


0.5 


n/2 


115 


0.323 


1.261 


0.753 


0.707 


0.188 




10^ 


1.0 


7T/2 


190 


0.203 


1.157 


0.649 


0.643 


0.184 


2- 


10^ 


0.5 


n/2 


131 


0.222 


1.162 


0.654 


0.643 


0.188 




10^ 


0.5 


n/2 


284 


0.306 


1.131 


0.623 


0.629 


0.177 




102 


0.5 


47T 


6 


0.703 


0.985 


0.478 


0.523 


0.201 


5- 


10^ 


0.5 


47T 


59 


0.237 


1.274 


0.766 


0.715 


0.188 




10^ 


0.5 


47T 


305 


0.159 


1.165 


0.657 


0.647 


0.184 



density. In general, the ignition process terminates after a few 
tenths of a second regardless of the exponentiation parameter 
(see Table [TJ. The case Cq = I corresponds to roughly one ig- 
nition per second and would fail to produce an explosion at all. 
On the other hand, there appears to be a remarkable robustness 
for Cq > 10^. The variation in the total nuclear energy release 
is less than 10% for all simulations with exponentiation pa- 
rameter in excess of 10^, while the number of ignition events 
may change by a factor as large as ten. Changing the ratio of 
the mixing length to the thermal radius has more or less the 
same effect as varying the exponentiation parameter. If R/A in- 
creases, the probability of ignition becomes higher. However, 
we found that choosing R/A significantly smaller than 0.5 in- 
troduces an artificial cutoff' in the range of ignition radii. For 
this reason, we set R/A = 0.5. 

The time evolution of the total energy is plotted for selected 
simulations in panel (a) of Fig. [l] One can clearly see the faster 
progression of the explosion with increasing number of igni- 
tions and the reduced rate of growth for Cq ~ 10^. The final 
kinetic energy of 7.5 • 10^^ erg in the case Cq ~ 10"^ is about 
75 % larger than for the reference simulation with a perturbed 
axisymmetric initial ffame (c3). The production of iron group 
and intermediate mass elements is shown in panels (b) and (d), 
respectively. It becomes clear that the statistical properties of 
the ignition process are crucial for the nucleosynthesis of iron 
group elements, whereas the the final amount of intermediate 
mass elements, which are represented by ^"^Mg, appears to be 
fairly robust. However, it is important to note that burning is 
terminated at the density threshold of lO^gcm"^ because we 
still lack a sensible treatment of the distributed burning regime 
during the late stage of the explosion. Continuing the burning 
process into this regime might entail significant changes in the 
predictions of intermediate mass element production. 

The integrated subgrid scale turbulence energy ^sgs is plot- 
ted in panel (c) of Fig. [I] In the early phase of the explosion, 
one can discern power-law growth of ^sgs, where the exponent 
appears to be universal for stochastic ignition. There is merely 
a temporal shift of the graph depending on the exponentiation 



parameter. From ^ 0.5 s onwards, the increase of SGS turbu- 
lence energy diminishes. This marks the onset of the fully tur- 
bulent regime. The peak of ^sgs is increasingly delayed as the 
ignition process progresses slower for smaller values of Cq. At 
still later time, the turbulence energy decreases due to the rapid 
bulk expansion and the quenching of the burning process which 
supplies energy to the ffow. The exceptionally high peak values 
of ^sgs for the less energetic explosions are caused by a major 
fraction of the released energy stirring the white dwarf matter 
into turbulent motion rather than being consumed by bulk ex- 
pansion. In this respect, the "dud explosions" obtained in the 
case of slow ingition are more like an overffowing boiling pot. 

With regard to the question whether deffagration models 
can account for observed type la supernovae, it is crucial to 
look at the distribution of nuclear species at a time when the 
explosion ejecta enter the phase of homologous expansion. To 
that end we performed three full-star simulations and analysed 
the ejecta at t = 5.0 s. From the statistics listed in Table [T] 
it might appear that the sensitivity to the ignition process is 
less pronounced than for the single-octant simulations, possibly 
due to the absence of artifical boundary conditions. However, 
Fig.|2H4|illustrate that there are, indeed, substantial diff'erences 
in the final outcome. For each simulation, the stratification of 
the mass density is shown in two-dimensional sections through 
the centre of the computational domain. It is instructive to com- 
pare the partial densities of unburned carbon and oxygen, iron 
group elements ("Ni") and intermediate mass elements ("Mg"). 
In the case Cq = 10^, the contour sections on the top of Fig.|2| 
show that a big lump of unburned material prevails in the cen- 
tral region. The corresponding fractional masses as function of 
radial velocity are plotted at the bottom of Fig. El Given the 
total amount of roughly 0.5 Mq of iron group elements, such 
an explosion definitely does not resemble any observed SN 
la. Increasing the exponentiation parameter, however, we ob- 
tain an ever more stratified distribution of nuclear species. For 
Cq = 5 • 10^, the nuclear ash appears to be mixed with unburned 
material over a wide range of radii and radial velocities, respec- 
tively (Fig. 13. In the case of most rapid ignition with Cq = 10^, 
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t [s] t [s] 



Fig. 1. Evolution of the total energy, the masses of iron group and intermediate mass elements, respectively, as well as the subgrid 
scale turbulence energy in simulations with axisymmetric initial flame (c3) and stochastic ignition with varying exponentiation 
parameter Ce, respectively. The computational domain is a single octant with N = 128^ grid cells. 



a clear trend towards concentration of nickel in the centre be- 
comes apparent, while carbon and oxygen remain mostly in 
the outer layers. Th is is difl'erent from the results reported by 
iRopke etall (l2QQ5l) for simultaneous ignitions, where the ma- 
terial at small radial velocities appears to be depleted of carbon 
and oxygen for a total number of initially burning bubbles of 
the order one hundred, while further increasing the number of 
bubbles reverses the trend. 

For a more detailed account of the flame evolution, we per- 
formed in addition one full- star simulation with the interme- 
diate exponentiation parameter Cq = 5-10^ using N = 512^ 
numerical cells, i.e. twice the initial spatial resolution. The ig- 
nition process for this simulation is illustrated in Fig.|5l It might 
seem unexpected that relatively few ignitions occur close to the 
centre. However, this is a natural consequence of the scaling of 
the ignition probability with the volume. Moreover, most of the 



material near the centre is burned quickly and so there remains 
little space for further ignitions. The bulk of ignition events 
takes place at radii in the range between 50 and 250 km from 
the centre. The relatively large radial spread might be an artifact 
of the statistical rather than local evaluation of the temperature 
excess ST\y in expression ^ for the inverse time scale of the 
Poisson process. In consequence, we definitely see the require- 
ment of an improved modelling of the temperature fluctuations. 
The conceptual design of a fluctuation model which allows for 
the computation of the spatiotemp oral distribution has recently 
been proposed bv lKersteinI (l2005l) . 

Three-dimensional snapshots of the flame fronts in the sim- 
ulation with high resolution are shown in Fig. |6l The subgrid 
scale turbulence velocity which determines the turbulent flame 
speed is indicated by contours at the flame surface. At time 
r = 0. 1 s (a), one can see a large number of small bubbles gener- 
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Fig. 2. Total and fractional mass densities in a central section Fig. 3. Total and fractional mass densities in a central section 
and the corresponding probability density functions in radial and the corresponding probability density functions in radial 
velocity space for a full-star simulation with Cq = 10^ at t = velocity space for a full-star simulation with Cq = 5 • 10^ at 



5.0 s. 



^ = 5.0s. 



ated by the stochastic ignition process. As the burning bubbles 
are floating outwards from the centre they begin to form the 
typical Rayleigh-Taylor mushroom shapes (b). At ^ = 0.5 s, the 
flames mostly have merged into a single structure which is sub- 
ject to intense turbulent flow (c). Apart from the small- scaling 
wrinkling of the flame front one can discern several dominat- 
ing modes. Particularly, a pronounced asphericity emerges after 
most of the explosive burning has come to pass (d). The forma- 
tion of dominating modes is a consequence of the fast growth 
of Rayleigh-Taylor instabilities seeded at large radii. Hence, the 
evolution of the flames is highly sensitive to the radial distribu- 
tion of the ignitions. This implication is particularly interesting 
in the light of recently found obse rvational evidence fo r as- 
phericity of type la supernova ejecta (iLeonard et al.l20o3) . The 
late-time evolution and, particularly, the distribution of nuclear 
species was found to be quite similar to the 256^ run. 



4. Conclusion 

We have achieved a first step toward the dynamical modelling 
of the ignition process in three-dimensional large-scale simula- 
tions of thermonuclear supernovae. There is one free parame- 
ter, called the exponentiation parameter Cq, which controls the 
statistical properties of the ignition process, especially, the to- 
tal number of ignition events. Varying Cq in the numerical sim- 
ulations alters the total mass of iron group elements and the 
stratification of burning products. Particular ignition events can 
aff'ect the morphology of the explosion ejecta as some modes 
grow faster than others. In fact, we found noticeable aspheric- 
ity of the developing flame front even for high values of Cq, 
where the statistical distribution of ignitions is nearly spheri- 
cally symmetric. 

If the total number of ignitions is less than one hundred in 
the whole core, only weak explosions with small total masses 
of iron group elements result. Although certain type la super- 
novae appear to produce, for instance, as little as 0.1 Mq of 
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Fig. 5. Plot of the radii of ignition events as function of time in 
the full- star simulation with N = 512^. 
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viously assumed. This would confirm the prediction of sim- 
ulations with an intermediate number of ignitions, say, about 
50 per octant (Fig.[3land|6|). Nevertheless, if the possibility of 
more than a few ignitions occurring within a tenth of a second 
was excluded by advances in understanding the pre- supernova 
core evolution, then the deflagration model would be ruled out 
as an explanation for typical observed SNe la. 
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Fig. 4. Total and fractional mass densities in a central section 
and the corresponding probability density functions in radial 
velocity space for a full-star simulation with Ce = 10^ at ^ = 
5.0 s. 



nickel (IStritzinger et al.lEo05l) . the exceptionally large amount 
of left-over carbon and oxygen found in the simulations with 
small number of ignitions renders them unlikely to correspond 
to real events. The picture might change, however, if burning 
at low density is properly taken into account and, possibly, 
substantially larger amounts of intermediate mass elements are 
produced. 

For a larger number of ignitions, on the other hand, strong 
explosions may result from pure deflagrations. However, even 
the most energetic explosions found in our simulations can 
not account for the upper range of observed nick el masses. In 
this re spect, we arrive at the same conclusion as iRopke etall 
(l2005l) for simultaneous multi-spot ignitions. On the other, 
we found that the explosion ejecta tend to exhibit an increas- 
ingly layered structure with growing number of ignition events. 
Consequently, the stochastic ignition model catches up with 
simulatio ns in which a delayed d etonation is initiated. Apart 
from that, 'Stritzinger et al.l (|2005") have found indications that 
possibly more nickel is mixed into the outer layers than pre- 
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(a) f = 0. 1 s, TVuni = 320^ Xuni = 815 km (b) ^ = 0.25 s, Nuni = 340^ X^ni = 1 125 km 




(c) t = 0.5 s, A^uni = 396^ Xuni = 2339 km (d) ^ = 0.75 s, A^uni = 438^ Xuni = 5177 km 



Fig. 6. Flame evolution in a full- star simulation with contours of the subgrid scale turbulence velocity in logarithmic scaling. 
For each snapshot, the physical time t and the number of cells, A^uni. as well as the size Xuni of the uniform part of the numerical grid which 

encloses the burning region are specified. 
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